
# Generate data


import numpy as np
import matplotlib.pyplot as plt
from math import pi, sin, cos, sqrt
from random import gauss, uniform, seed
import sys
import csv


Gkp = 0
Gki = 0
Gkd = 0



#=============================================================================
Gintegral = 0
Ge0 = 0

#-----------------------------------------------------------
# Sub-rotina de controle.
def pid (e1 \
       ,kp = 0.00001 \
       ,ki = 0.0001 \
       ,kd = 0.005) :

   ctrl = (e1        * kp) +  \
         (Gintegral  * ki) +  \
         ((e1 - Ge0) * kd)

   Gintegral = Gintegral + 0.5 * (Ge0 + e1)

   Ge0 = e1

   return ctrl




#=============================================================================

# For reproducibility purposes
rseed = 123
seed(rseed)


GReagenteTotal    = 0.0
GConcentraEntrada = 1.0      # gramas por litro
GVolumeEntrada    = 10.0     # litros
GVolumeSaida      = 10.0     # litros
GVolumeReator     = 1000.0   # litros
GTempo            = 0


#-----------------------------------------------------------
def Valve(t):

   return gauss(1.0, 0.01)  * \
          1.0 + ((1.25 ** (t/100.0)) / 5000.0)


#-----------------------------------------------------------
def UpdateReator():

   global GReagenteTotal
   global GConcentraEntrada
   global GVolumeEntrada
   global GTempo

   GReagenteTotal = GReagenteTotal + \
                    Valve(GTempo) * GVolumeEntrada * GConcentraEntrada - \
                    GVolumeSaida  * (GReagenteTotal / GVolumeReator) 

 
#-----------------------------------------------------------


nivel = []

for GTempo in range(0,2000):

   UpdateReator ()

   nivel.append (GReagenteTotal)


try:
   with open("Reference.csv", "w", newline='') as arqout : 
      escritorcsv    = csv.writer(arqout)
      escritorcsv.writerow (nivel)

   print ("Trajetória de referência gravada em <Reference.csv>.")

except:
  print ("Erro ao gravar dados medidos.")

quit()
  

plt.plot (nivel)
plt.title  ('Reference', fontsize=14)
plt.xlabel ('Time', fontsize=14)
plt.ylabel ('Solute concentration', fontsize=14)

ax = plt.gca()

#ax.set_xlim([xmin, xmax])
ax.set_ylim([950, 1050])

for x in (\
          #[ax.title,ax.xaxis.label, ax.yaxis.label] + \
          ax.get_xticklabels() + ax.get_yticklabels()):
    x.set_fontsize(12)

plt.grid ()
plt.show ()


